Climate Validation

1 Packages

library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(Hmisc) #for correlation matrix

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

2 Read Data

NA.month <- read.csv("ClimateNA/outputs/monthly_tidy.csv")
inst.month.P <- read.csv("Instrumental/outputs/monthly_precipitation.csv")
inst.month.T <- read.csv("Instrumental/outputs/monthly_temperature.csv")

3 Prepare Data

NA.month$Month <- as.factor(NA.month$Month)
inst.month.P$month <- as.factor(inst.month.P$month)
inst.month.T$month <- as.factor(inst.month.T$month)

m.m.t <- NA.month %>% select(Year, Month, Tave)

i.m.t <- inst.month.T %>% 
  select(year, month, mean_temp_year) %>%
  filter(year > 1900) %>%
  rename(Year = year,
         Month = month,
         Tave = mean_temp_year)

m.m.t$Source <- "ClimateNA"
i.m.t$Source <- "Instrumental"

months.T <- rbind(m.m.t, i.m.t)

m.m.p <- NA.month %>% select(Year, Month, PPT)

i.m.p <- inst.month.P %>% 
  select(year, month, total_precip_year) %>%
  filter(year > 1900) %>%
  rename(Year = year,
         Month = month,
         PPT = total_precip_year)

m.m.p$Source <- "ClimateNA"
i.m.p$Source <- "Instrumental"

months.P <- rbind(m.m.p, i.m.p)

4 Visualize Data

4.1 Months

ggplot() + 
  geom_boxplot(data = NA.month, aes(x = Month, y = Tave), alpha = 0.2, col = "darkblue", notch = TRUE) + 
  geom_boxplot(data = inst.month.T, aes(x = month, y = mean_temp_year), alpha = 0.2, col = "darkgreen", notch = TRUE) + 
  theme_bw() +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

ggplot() + 
  geom_boxplot(data = NA.month, aes(x = Month, y = PPT), alpha = 0.2, col = "darkblue", notch = TRUE) + 
  geom_boxplot(data = inst.month.P, aes(x = month, y = total_precip_year), alpha = 0.2, col = "darkgreen", notch = TRUE) + 
  theme_bw() +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

ggplot(data = months.T, aes(x=Year, y = Tave, col = Source)) +
  geom_line() +
  facet_grid(Month~.) + 
  theme_classic() 

ggplot(data = months.P, aes(x=Year, y = PPT, col = Source)) +
  geom_line() +
  facet_grid(Month~.) + 
  theme_classic() 

# ggplot() +
#  geom_line(data = study_p, aes(x = month, y = total_precip), col= "skyblue", size = 1) +
#   geom_ribbon(data = study_p, aes(x = month, ymin = (lower), ymax = (upper)), fill = "skyblue", alpha = 0.2) +
#   geom_line(data = NA_p, aes(x = Month, y = mean_precip), col= "steelblue4", size = 1) +
#   geom_ribbon(data = NA_p, aes(x = Month, ymin = (lower), ymax = (upper)), alpha = 0.2, fill= "steelblue4") +
#   scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10,11,12)) +
#   theme_bw() +
#    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

4.2 Annual

# ggplot() +
#   # London Data
#   geom_line(data=london_temp_annual, aes(x = year, y = mean_temp), colour = "lightseagreen", linetype = "dashed", alpha = 0.8) +
#   # Pinery Data
#   geom_line(data=pinery_temp_annual, aes(x = year, y = mean_temp), colour="lightsalmon4", linetype = "dashed", alpha = 0.8) +
#   # ClimateNA
#   geom_line(data=annual, aes(x = Year, y = MAT), colour="darkgreen") +
#   #Aesthetics
#   xlab("Year") +
#   ylab("Mean Annual Temperature (°C)") +
#   scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
#   #ggtitle("Average Annual Temperature in \n Southwestern Ontario from 1883 to 2017") + 
#   theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
#  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

5 Model

i.m.p <- i.m.p %>% rename(Instrumental=PPT) %>% select(-Source)
m.m.p <- m.m.p %>% rename(ClimateNA=PPT)  %>% select(-Source)

p.complete <- right_join(i.m.p, m.m.p, by = c("Year", "Month")) %>% na.omit()
  
p.matrix <- p.complete %>% select(-Year, -Month)

i.m.t <- i.m.t %>% rename(Instrumental=Tave) %>% select(-Source)
m.m.t <- m.m.t %>% rename(ClimateNA=Tave)  %>% select(-Source)

t.complete  <- right_join(i.m.t, m.m.t, by = c("Year", "Month")) %>% na.omit()

t.matrix <- t.complete %>% select(-Year, -Month)

5.1 Fit Model

m.t1 <- lm(data = t.complete, ClimateNA ~ Instrumental)
summary(m.t1)

## 
## Call:
## lm(formula = ClimateNA ~ Instrumental, data = t.complete)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.41302 -0.35540  0.00878  0.33167  1.85683 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.418485   0.018835   22.22   <2e-16 ***
## Instrumental 0.978125   0.001511  647.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5205 on 1246 degrees of freedom
## Multiple R-squared:  0.997,  Adjusted R-squared:  0.997 
## F-statistic: 4.192e+05 on 1 and 1246 DF,  p-value: < 2.2e-16

Remove outliers

p.complete <- p.complete %>% filter(Instrumental < 200)

Relevel

p.complete$Month <- relevel(p.complete$Month, ref = 5)

m.p1 <- lm(data = p.complete, ClimateNA ~ Instrumental)
summary(m.p1)

## 
## Call:
## lm(formula = ClimateNA ~ Instrumental, data = p.complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -54.940 -11.719  -0.927  10.306  76.853 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   27.1542     1.2108   22.43   <2e-16 ***
## Instrumental   0.5411     0.0141   38.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.46 on 1234 degrees of freedom
## Multiple R-squared:  0.5439, Adjusted R-squared:  0.5435 
## F-statistic:  1472 on 1 and 1234 DF,  p-value: < 2.2e-16

m.p2 <- lm(data = p.complete, ClimateNA ~ Instrumental*Month)
summary(m.p2)

## 
## Call:
## lm(formula = ClimateNA ~ Instrumental * Month, data = p.complete)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.608 -11.389  -1.025  10.188  79.701 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          19.680204   3.886351   5.064 4.74e-07 ***
## Instrumental          0.655126   0.043422  15.087  < 2e-16 ***
## Month1               10.222204   5.888583   1.736 0.082829 .  
## Month2                2.832610   5.844245   0.485 0.627988    
## Month3                2.557246   5.862190   0.436 0.662750    
## Month4                5.253099   6.141415   0.855 0.392523    
## Month6               11.500457   5.634142   2.041 0.041446 *  
## Month7               21.570780   5.249084   4.109 4.23e-05 ***
## Month8               11.509542   5.343394   2.154 0.031439 *  
## Month9                1.337119   5.507942   0.243 0.808231    
## Month10              -1.132427   5.649630  -0.200 0.841168    
## Month11              -0.894239   6.150210  -0.145 0.884419    
## Month12              15.112783   6.166552   2.451 0.014396 *  
## Instrumental:Month1  -0.180905   0.066507  -2.720 0.006620 ** 
## Instrumental:Month2  -0.096321   0.076943  -1.252 0.210869    
## Instrumental:Month3  -0.062959   0.076407  -0.824 0.410105    
## Instrumental:Month4   0.006738   0.073945   0.091 0.927413    
## Instrumental:Month6  -0.081111   0.065660  -1.235 0.216955    
## Instrumental:Month7  -0.262228   0.058347  -4.494 7.65e-06 ***
## Instrumental:Month8  -0.204063   0.059557  -3.426 0.000632 ***
## Instrumental:Month9  -0.029312   0.060796  -0.482 0.629792    
## Instrumental:Month10 -0.017195   0.064328  -0.267 0.789279    
## Instrumental:Month11 -0.021626   0.066670  -0.324 0.745708    
## Instrumental:Month12 -0.266745   0.066686  -4.000 6.72e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.84 on 1212 degrees of freedom
## Multiple R-squared:  0.5832, Adjusted R-squared:  0.5753 
## F-statistic: 73.73 on 23 and 1212 DF,  p-value: < 2.2e-16

t.complete$fit.1 <- fitted(m.t1) #Model predictions or fitted values
p.complete$fit.1 <- fitted(m.p1) 
p.complete$fit.2 <- fitted(m.p2)

5.2 Model Predictions

ggplot(data = t.complete) +
  geom_line(aes(x=Instrumental, y = fit.1, col = Month))  +
  geom_point(aes(x=Instrumental, y = ClimateNA, col = Month), size = 0.5, alpha = 0.3) +
     labs(x = "Instrumental: Average Monthly Temperature",
       y = "ClimateNA: Average Monthly Temperature") +
  theme_classic()

Relevel to plot

p.complete$Month <- factor(p.complete$Month, levels = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12"))

ggplot(data = p.complete) +
  geom_line(aes(x=Instrumental, y = fit.1, col = Month))  +
  geom_point(aes(x=Instrumental, y = ClimateNA, col = Month), size = 0.5, alpha = 0.5) +
   labs(x = "Instrumental: Monthly Precipitation (mm)",
       y = "ClimateNA: Monthly Precipitation (mm)") +
  theme_classic()

ggplot(data = p.complete) +
  geom_line(aes(x=Instrumental, y = fit.2, col = Month))  +
  geom_point(aes(x=Instrumental, y = ClimateNA, col = Month), size = 0.5, alpha = 0.4) +
  lims(x = c(0,200), y = c(0,200)) +
   labs(x = "Instrumental: Monthly Precipitation (mm)",
       y = "ClimateNA: Monthly Precipitation (mm)") +
  theme_classic()

5.3 Assumptions

Normality:

plot(m.t1, which=2)

plot(m.p1, which=2)

plot(m.p2, which=2)

Equal variance:

plot(m.t1, which=3)

plot(m.p1, which=3)

plot(m.p2, which=3)

Linearity:

plot(m.t1, which=1)

plot(m.p1, which=1)

plot(m.p2, which=1)

Leverage:

plot(m.t1, which=4)

plot(m.p1, which=4)

plot(m.p2, which=4)

6 Save Outputs

7 Reproducibility

Sys.time()

## [1] "2020-06-22 16:32:42 PDT"

git2r::repository()

## Local:    master /Users/JenBaron/Documents/UWO/UWO NSERC/Statistical Analysis/Climate/Climate_Pinery
## Remote:   master @ origin (https://github.com/JenBaron/Climate_Pinery.git)
## Head:     [ad4302f] 2020-06-19: seasonal climateNA

sessionInfo()

## R version 4.0.0 (2020-04-24)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.5
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] Hmisc_4.4-0     Formula_1.2-3   survival_3.1-12 lattice_0.20-41
## [5] gridExtra_2.3   tidyr_1.0.2     dplyr_1.0.0     ggplot2_3.3.1  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0    xfun_0.14           purrr_0.3.4        
##  [4] splines_4.0.0       colorspace_1.4-1    vctrs_0.3.1        
##  [7] generics_0.0.2      htmltools_0.4.0     yaml_2.2.1         
## [10] base64enc_0.1-3     rlang_0.4.6         pillar_1.4.4       
## [13] foreign_0.8-79      glue_1.4.1          withr_2.2.0        
## [16] RColorBrewer_1.1-2  jpeg_0.1-8.1        lifecycle_0.2.0    
## [19] stringr_1.4.0       munsell_0.5.0       gtable_0.3.0       
## [22] htmlwidgets_1.5.1   evaluate_0.14       labeling_0.3       
## [25] latticeExtra_0.6-29 knitr_1.28          htmlTable_1.13.3   
## [28] Rcpp_1.0.4.6        acepack_1.4.1       scales_1.1.1       
## [31] backports_1.1.7     checkmate_2.0.0     farver_2.0.3       
## [34] png_0.1-7           digest_0.6.25       stringi_1.4.6      
## [37] grid_4.0.0          tools_4.0.0         magrittr_1.5       
## [40] tibble_3.0.1        cluster_2.1.0       crayon_1.3.4       
## [43] pkgconfig_2.0.3     ellipsis_0.3.1      Matrix_1.2-18      
## [46] data.table_1.12.8   rmarkdown_2.1       rstudioapi_0.11    
## [49] R6_2.4.1            rpart_4.1-15        git2r_0.26.1       
## [52] nnet_7.3-14         compiler_4.0.0